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Different relaxation approximations to partial differential equations, including 
conservation laws, Hamilton- Jacobi equations, convection-diffusion problems, 
gas dynamics problems, have been recently proposed. The present paper focuses 
onto diffusive relaxed schemes for the numerical approximation of nonlinear re- 
action diffusion equations. High order methods are obtained by coupling ENO 
and WENO schemes for space discretization with fMEX schemes for time in- 
tegration, where the implicit part can be explicitly solved at a linear cost. To 
illustrate the high accuracy and good properties of the proposed numerical 
schemes, also in the degenerate case, we consider various examples in one and 
two dimensions: the Fisher-Kolmogoroff equation, the porous-Fisher equation 
and the porous medium equation with strong absorption. 
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1. Introduction 

The main purpose of this work is to approximate solutions of a nonlinear, 
possibly degenerate, reaction-diffusion equation of the form 

^ = DA(p(u))+g(u) (1) 

for x £ f2 C R d , d > 1, t > 0, with initial condition, u(x, 0) = Uq{x) and 
with suitable boundary conditions, where the function p is a non-decreasing 
Lipschitz continuous function defined on R. The equation is degenerate if 
p(0) = and in this case the solutions often become non-smooth in finite 
time, developing fronts and discontinuities [1]. Finally, the coefficient D is 
a diffusivity coefficient and the function g(u) is the reaction term. 

Equations like (1) are relevant in describing biological and physical pro- 
cesses and are involved in the modelling of population growth and dispersal, 
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of waves of concentration of chemical substances in living organisms, of the 
motion of viscous fluids (see e.g. [2]). Also, systems of equations like (1), 
coupled via the reaction terms, are capable of modelling the cyclic Belousov- 
Zhabotinskii reactions and the pattern formation on the wings of butterflies 
and the coat of mammals (see [3] and references therein). 

This paper is organized as follows. In section 2 we introduce the relax- 
ation approximation of nonlinear diffusion problems, in section 3 we de- 
scribe the fully discrete relaxed numerical scheme in the reaction-diffusion 
case. In section 4 we report several numerical tests, both in one and two 
space dimensions. 

2. Relaxation approximation of nonlinear diffusion 

The schemes proposed in the present work are based on the same frame- 
work of the well-known relaxation approximation of [4] for hyperbolic con- 
servation laws. In the case of the nonlinear diffusion operator in (1), an 
additional variable v(x, t) G R d and a positive parameter e are introduced 
and the following relaxation system is obtained: 



Now, formally, in the small relaxation limit, e — > + , system (2) approxi- 
mates to leading order equation (1). The parameter ip is introduced in order 
to have bounded characteristic velocities and to avoid a singular differential 
operator as e — > + . 

Finally the non linearity in the convective term is removed, as in stan- 
dard relaxation schemes, introducing another scalar variable w(x,t) and 
rewriting the system as: 



Formally, as e — > + , w — > p(u), v — ► — Vp(u) and the original equation is 
recovered. 



— + div(u) = g(u) 




(2) 




(3) 
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In the previous system the parameter e has physical dimensions of time 
and represents the so-called relaxation time. Furthermore, w has the same 
dimensions as u, while each component of v has the dimension of u times 
a velocity; finally ip is a velocity. The inverse of e gives the rate at which v 
decays onto — Vp(u) in the evolution of the variable v governed by the stiff 
second equation of (3). 

Equations (3), originally introduced in [5] for the purely diffusive case, 
form a semilinear hyperbolic system with a stiff source term. The charac- 
teristic velocities of the hyperbolic part are given by 0, ±<p. The parameter 
ip allows to move the stiff terms ®Wp{u) to the right hand side, without 
losing the hyperbolicity of the system. 

We point out that degenerate parabolic equations often model physical 
situations where free boundaries and discontinuities are relevant: we expect 
that schemes for hyperbolic systems will be able to reproduce faithfully 
these details of the solution. One of the main properties of (3) consists in the 
semilinearity of the system, that is all the nonlinear terms are in the (stiff) 
source terms, while the differential operator is linear. Hence, the solution of 
the convective part requires neither Riemann solvers nor the computation 
of the characteristic structure at each time step, since the eigenstructure 
of the system is constant in time. Moreover, the relaxation approximation 
does not exploit the form of the nonlinear function p and hence it gives rise 
to a numerical scheme that, to a large extent, is independent of it, resulting 
in a very versatile tool. 

We also anticipate here that, in the relaxed case (i.e. e = 0), the stiff 
source terms can be integrated solving a system that is already in triangular 
form and then it does not require iterative solvers. 



3. Relaxed numerical schemes 
3.1. Relaxed IMEX schemes 

We extend here the schemes studied in [6] , including the reaction term. For 
simplicity, here we describe the one-dimensional case, the generalization 
being straightforward. 

We observe that system (3) is in the form 

dz df(z) , . 1 , , , 

where z = (u,v,w) T 7 f(z) = (v,ip 2 w,v) T , g(z) — (g(u), 0, 0) T and h(z) — 
(0, — v + (e(p 2 — D)w x ,p(u) — w) T . When e is small, the presence of both 
non-stiff and stiff terms, suggests the use of IMEX schemes [7-9]. In the 
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problems considered here, the reaction term g(u) is not stiff and hence we 
treat it with the explicit portion of the IMEX scheme. 

First we consider a semi-discrete form for time integration. Let's assume 
for simplicity a uniform time step At and denote with z n the numerical 
approximation of the variable z at time t n = nAt, for n = 0, 1, ... In our 
case a ^-stages IMEX scheme reads 



z n+1 = z n - AtJ2~b t 



i=l 



dx 



1=1 



where the stage values are computed as 

At 



z« =sW+ a iti h(z^) 



for 



i-i 



fc=i 



^(^)+ 9 (zW) 



dx 



+ ^ i £a iik h(zW) 
£ k=i 



(4) 



(5) 



(6) 



Here (aik,bi) and (dik,bi) are a pair of Butcher's tableaux [10] of, respec- 
tively, a diagonally implicit and an explicit Runge-Kutta schemes. 

In this work we use the so-called relaxed schemes, that are obtained by 
formally letting e — > in (4). For these the computation of the first stage, 
that is (5) with i = 1, 



w 



(1) 



u" 
w % 



— ai \h 

£ 



W (D 



(7) 



implies that h(z^) = 0, which is equivalent to 



dx 



Now the second stage, i = 2, reads 
^) = ,"-Ate 2 , 1 fg(z( 1 ))+.g(zW) 



+ _^ a21 / l(z (i)) + _^ a22 / l(z (2) ) . 



Hence the last two components of z^ are determined by the stiff terms of 
the above expression, namely h(z^) — 0. On the other hand, due the form 
of h(z), there are no stiff terms in the equation for the first component , 
which is then determined by a balance law. 
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Summarizing, the relaxed scheme yields an alternation of relaxation 
steps 

h(z^) = i.e. w w = v« = 

and transport steps where we advance for time a^Ai 

with initial data z — z^ l \ retain only the first component and assign it to 
pi na Uy the value of u n+1 is computed as u n + hu^K 



3.2. Spatial reconstructions 

In order to have a fully discrete scheme, we still need to specify the space 
discretization. We will use discretizations based on finite differences, in 
order to avoid cell coupling due to the source terms. 

Recall that the IMEX technique reduces the integration to a cascade of 
relaxation and transport steps. The former are the implicit parts of (7) and 
(5), while the transport steps appear in the evaluation of the explicit terms 
£?W in (6). Since (7) and (5) involve only local operations, the main task 
of the space discretization is the evaluation of d x f, where we will exploit 
the linearity of / in its arguments. 

In the one-dimensional case, let us consider a uniform grid on [a, b] C K, 
xj = a — -| + jh for j = 1, . . . , m, where h — (b — a)/m is the grid spacing 
and m the number of cells. We denote with z 3 n the value of the quantity z 
at time t n at Xj , the centre of the j th computational cell. The fully discrete 
scheme may be written as 

- - A* ± h (F?l /2 - Ffl 1/2 ) + ^ ± b i9 (zf), 
1=1 1=1 

where F^ 1 }. ,„ are the numerical fluxes, which are the only item that we still 
need to specify. It is necessary to write the scheme in conservation form 
and thus, following [11], we introduce the function F such that 

-, ,-x+h/2 

f(z(x,t)) = - F(s,t)ds 

11 Jx-h/2 

and hence 

^(z( Xj ,t)) = i (P{x J+1/2 ,t) - F( Xj _ 1/2 ,tfj . 
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The numerical flux function approximates F(x J+1 / 2 )- 

In order to compute the numerical fluxes, for each stage value, we re- 
construct boundary extrapolated data z j+i/ 2 with a non-oscillatory inter- 
polation method, starting from the point values of the variables at 
the centre of the cells. Next we apply a monotone numerical flux to these 
boundary extrapolated data. 

To minimize numerical viscosity we choose the Godunov flux, which, 
in the present case of a linear system of equations, reduces to the upwind 
flux. In order to select the upwind direction we write the linear system with 
constant coefficients (8) in characteristic form. The characteristic variables 
relative to the eigenvalues ip,—ip,0 are respectively 

v + tpw tpw-v 
U = V = W = u — w. 

Note that u = U + V + W . Therefore the numerical flux in characteristic 
variables is F j+1/2 = {<pUj~ +1/2 ,-<pV£_ 1/2 ,Q). 

The accuracy of the scheme depends on the accuracy of the recon- 
struction of the boundary extrapolated data. For a first order scheme 
we use a piecewise constant reconstruction such that Uj~ +1 i 2 = Uj and 
Vj^ 1 / 2 = Vj+i- F° r higher order schemes, we use ENO or WENO recon- 
structions of appropriate accuracy [11]. 

Since the transport steps need to be applied only to uW, we have 

i-l 

fe=i 

Finally, taking the last stage value and going back to conservative variables, 

u] +1 = «? - 1 Er=i bi 027/2 + «jjt /a - {vf- /3 + vf_x /2) ] 

+^ /2 - W ^ /2 -(w^ /2 - W ^ /2 )\) 

We wish to emphasize that the scheme reduces to the time advancement 
of the single variable u. Although the scheme is based on a system of three 
equations, the construction is used only to select the correct upwinding for 
the fluxes of the relaxed scheme and the computational cost of each time 
step is not affected. 

3.3. Numerical scheme 

Employing a Runge-Kutta IMEX scheme of order p, can give an integration 
procedure for (1) which is of order up to 2p with respect to h (see [6]), 



wiU {k) ~ -U {k) ~ -V {k)+ + V [k)+ 1 + Ata(u (k) ) 

^1^+1/2 U j-l/2 V j+l/2 + V j-l/2 J + H9K U j ) 
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because the CFL restrictions on the time-step are of parabolic type (At < 
Ch 2 ). We observed that this theoretical convergence rate can be achieved in 
practice, with careful choice of the approximations of the spatial derivatives. 

Summarizing, the relaxed schemes that we propose for the numerical 
integration of (1) consists of the following steps. 

For each Runge-Kutta stage we need to compute the variables and 
tyW (relaxation steps). The computation of t>W requires the approximation 
of a spatial gradient operator, for which we choose a central finite difference 
operator of order at least 2p. Then we need to solve the transport equation 
by diagonalizing the linear system (8), reconstructing the characteristic 
variables 17 W and at cell boundaries and computing the fluxes. Again, 
we must choose a spatial reconstructions of order at least 2p and, to avoid 
spurious oscillations, we employ ENO or WENO non-oscillatory procedures. 
These procedures compare, for each cell, the reconstructions obtained using 
different stencils and choose the least oscillatory one (ENO) or compute a 
weighted linear combination of all of them (WENO). For details, see [11]. 

4. Numerical results 

4.1. Travelling waves tests 

As a first test, we consider the one-dimensional Fisher-Kolmogoroff equa- 
tion, namely 



for x € K and t > 0, which was initially proposed for the modelling of 
spatial spread of populations [2] . 

The two uniform steady states of (9) are uo(x) = and u\(x) = 1. A 
careful analysis of the state space, reveals the existence of travelling wave 
solutions linking these two states, i.e. of solutions of the form u(x, t) = 
U(x — ct) such that U{— oo) = 1 and U(oo) — 0. [2] reports the following 
asymptotic form for such solutions: 



which is valid when the speed verifies 

c > 2VkD. 

Moreover the slope of the inflection point (z = 0) is related to the speed c 




(9) 




(10) 
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c=2t=0,2,4,6,8,10 




-SO SO 100 150 200 



Fig. 1. Travelling waves with c = 2 (top) and c = 10 (bottom). The crosses represent 
the numerical solution at t = 0, 2, 4, 6, 8, 10 and the solid lines are the corresponding 
asymptotic expansions (10). 

of the wave by the relation 

C/'(0) = -l + o(l),forc>2. (11) 

Travelling wave solutions are very important in the modelling of popula- 
tions, since they represents phenomena like the invasion of a territory by a 
new species, the expansion of epidemics, etc. 

We tested the ability of our integration scheme to reproduce these re- 
sults. We chose k = D = 1 in (9), set up initial datum u(x, 0) = U(x) 
for c = 2, 4, . . . , 10 on x € [—50,200] and evolved it with homogeneous 
Neumann boundary conditions until T = 10. We present the numerical so- 
lutions obtained with ENO spatial reconstruction of order 6 and third order 
Runge-Kutta time integration: they are compared with the asymptotic ones 
in Figure 1. Moreover we tested the validity of the expansion (11) plotting 
the maximum of the numerical gradient of the solution, together with the 
values predicted by the asymptotic analysis. From Figure 2 we clearly see a 
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Fig. 2. Comparison of the maximum slope of the numerical solutions (dots) with the 
value predicted by the asymptotic expansion (11). Time is on the horizontal axis, the 
values of c are printed on the graph close to each line. 



very accurate agreement, except for the limit case c = 2, which is however 
at the boundary of the validity range of (11). 



1=0.00, 0.30, 0.62, 0.92, 1.24, 1 .55, 1.87, 2.27, 2.65, 3.05, 3.43, 3.83 




-8 -G -4 -2 2 4 6 



Fig. 3. Comparison of the numerical solutions (dots) of (12) with the asymptotic ex- 
pansion (solid line). Here p = 1 and q = m. The latter is printed only in the time range 
where it is valid. 

As a second test we consider the following generalization of the Fishcr- 
Kolmogoroff equation (9) 

(--I) < 12 » 

The existence of travelling waves can be proved for a wide range of param- 
eters p,q,m [2]. The paper [12] gives an expression of such waves for the 
case p = 1 and q — m. Moreover [12] finds an asymptotic expansion for the 
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case of two merging travelling waves, which is valid for a finite time after 
the first contact of the two fronts. 

Figure 3 shows both the numerical solutions obtained with our method 
and the asymptotic expansion given in [12], in the case q = m = 1. The 
initial data represent two travelling fronts, initially at xq = — 1 and x\ = 1, 
moving in opposite directions. They first meet at t* = 1.41 and then the 
two waves merge. The asymptotic expansion given in [12] is valid for a 
small time interval after t* and it is shown in Figure 3 with solid lines. The 
numerical solutions (in dots) are shown until the steady state is reached. 

4.2. Two dimensional tests 



1=000 1=007 




Fig. 4. The numerical solution of (13) with initial data (14) at different times, for 
(x,y) G [-1.2, 1.2] 2 . 



We tested the multi-dimensional version of our integration scheme on 
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the equation 



du 
~dt 



= A(u m ) - cvP 



(13) 



for x G [—2, 2] 2 and t > 0. The above equation presents interesting finite- 
time extinction phenomena, reported in [13]. We first take initial data 



1 



u(x,y,0) 



1 - 



\/x 6 +y 6 



,x = 0,y = 
otherwise 



(14) 



and evolve it until extinction with equation (13) for c = 5, p = 0.5 and 
m = 2. The results are presented in Figure 4. 

Finally we tested the persistence of asymmetry in the initial datum along 
the evolution. We took u(x,y,0) as a radially symmetric function with a 
small perturbation, see Figure 5, and evolved it with the same equation 
and parameters as before, until extinction. Figure 5 shows clearly that the 
initial perturbation of the radial symmetry is maintained until the solution 
vanishes. 
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5. Conclusions 

We have proposed and analyzed relaxed schemes for nonlinear degenerate 
reaction diffusion equations. By using suitable discretization in space and 
time, namely ENO/WENO non-oscillatory reconstructions for numerical 
fluxes and 1MEX Runge-Kutta schemes for time integration, we have ob- 
tained a class of high order schemes. The theoretical convergence analysis 
for the semidiscrete scheme and the stability for the fully discrete schemes 
have been studied by us, for the case of nonlinear diffusion, in [6] . 

Here we tested these schemes on travelling waves solutions and in cases 
where the solution vanishes in finite time. In all cases we observed a very 
good agreement with known properties of the exact solutions. 
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